##########################################################################################

library(data.table)
library(optparse)
library(Seurat)
library(ArchR)
library(ggthemes)
library(dplyr)
library(ggsci)

##########################################################################################
option_list <- list(
    make_option(c("--rna_file"), type = "character"),
    make_option(c("--comine_data_all_file"), type = "character"),
    make_option(c("--gene_list_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    ## 单细胞表达文件
    rna_file <- "~/20231121_singleMuti/results/qc_atac_v3/germ/testis_combined.annotationCellType.qc.Rdata"

    ## 所有的细胞的,计算maxdelt
    comine_data_all_file <- "~/20231121_singleMuti/results/qc_atac_v3/germ/testis_combined_peak.combineRNA.qc.Rdata"

    ## 认为关键的TF
    gene_list_file <- "~/20231121_singleMuti/config/Human_reported_TF2.csv"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/report_tf_germ"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

rna_file <- opt$rna_file
comine_data_all_file <- opt$comine_data_all_file
gene_list_file <- opt$gene_list_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################

grp_order2 = c("SSC",
"Differenting&Differented SPG",
"Leptotene",
"Zygotene",
"Patchytene",
"Diplotene",
"Early stage of spermatids",
"Round&ElongateS.tids",
"Sperm",
"Leydig cells",
"Myoid cells",
"Pericytes",
"Sertoli cells",
"Endothelial cells",
"NKT cells",
"Macrophages")

use_colors <- c(pal_npg("nrc")(10) , pal_jco("default")(6))
#names(use_colors) <- unique(scrnat$cell_type)
names(use_colors) <-c( "Myoid cells","Leydig cells" ,          
"Endothelial cells","Zygotene",         
"Round&ElongateS.tids","Patchytene",
"SSC","Sperm" ,                 
"Diplotene","Early stage of spermatids", 
"Leptotene","Sertoli cells",            
"Macrophages","Differenting&Differented SPG",
"Pericytes","NKT cells")

names(use_colors) <- paste0( as.numeric(factor( names(use_colors) , levels = grp_order2 , order = T )) , names(use_colors))

###########################################################################################
## 导入数据
a <- load(rna_file)
## scrnat
DefaultAssay(scrnat) <- "MAGIC_RNA"
scrnat$cell_type <- factor( scrnat$cell_type , levels = grp_order2[grp_order2 %in% unique(scrnat$cell_type)] , order = T )
## 确认排序
Idents(scrnat) <- scrnat$cell_type

b <- load(comine_data_all_file)
## testis_combined_peak_combineRNA
projHeme5 <- testis_combined_peak_combineRNA

## 细胞类型排序
projHeme5@cellColData$cell_type <- factor( projHeme5@cellColData$cell_type , levels = grp_order2 , order = T )
projHeme5@cellColData$cell_type_plot <- paste0(as.numeric(projHeme5@cellColData$cell_type) , projHeme5@cellColData$cell_type)

gene_list <- data.frame( fread(gene_list_file , header = F) )$V1

###########################################################################################
## rna测到的基因
rna_gene <- rownames(scrnat@assays$MAGIC_RNA)
## atac存在开放程度的基因
atac_gene <- getFeatures(projHeme5, useMatrix = "GeneScoreMatrix")
exp_gene <- getFeatures(projHeme5, useMatrix = "GeneExpressionMatrix")

##得到motif的矩阵
motif_gene <- getFeatures(projHeme5 , useMatrix = "MotifMatrix")
motif_gene <- sapply(strsplit( sapply(strsplit(grep( "z:" , motif_gene , value = T ) , ":") , "[" , 2) , "_" ) , "[" , 1)

###########################################################################################

#gene_list <- gene_list[gene_list[gene_list[gene_list %in% rna_gene] %in% atac_gene] %in% motif_gene]
gene_list <- gene_list[gene_list %in% c(rna_gene , atac_gene , motif_gene , exp_gene)]

## 如果三者均存在
for( gene in gene_list ){

    if(gene %in% motif_gene){
        markergene <- getFeatures(projHeme5, select = paste(gene, collapse="|"), useMatrix = "MotifMatrix")
        markergene <- grep("z:", markergene, value = TRUE)

        ## motif的活性
        pointSize <- 1
        p1 <- plotEmbedding(
            ArchRProj = projHeme5, 
            colorBy = "MotifMatrix", 
            name = sort(markergene), 
            embedding = "UMAP",
            imputeWeights = getImputeWeights(projHeme5),
            plotAs="points", size = pointSize
        )
        p1_violin <- plotGroups(
            ArchRProj = projHeme5, 
            groupBy = "cell_type_plot",
            colorBy = "MotifMatrix", 
            name = sort(markergene), 
            imputeWeights = getImputeWeights(projHeme5),
            plotAs="violin", size = pointSize ,
            useGroups = grp_order2,
            pal = use_colors ,
            alpha = 0.4,
            log2Norm = TRUE,
            addBoxPlot = TRUE
        )
    }

    if(gene %in% atac_gene){
        ## TF的活性
        p2 <- plotEmbedding(
            ArchRProj = projHeme5, 
            colorBy = "GeneScoreMatrix", 
            name = sort(gene), 
            embedding = "UMAP",
            quantCut = c(0.01, 0.95),
            imputeWeights = getImputeWeights(projHeme5),
            pal = paletteContinuous(set = "solarExtra")
        )

        p2_violin <- plotGroups(
            ArchRProj = projHeme5, 
            groupBy = "cell_type_plot",
            colorBy = "GeneScoreMatrix", 
            name = sort(gene), 
            imputeWeights = getImputeWeights(projHeme5),
            plotAs="violin", size = pointSize ,
            pal = use_colors ,
            useGroups = grp_order2,
            alpha = 0.4,
            log2Norm = TRUE,
            addBoxPlot = TRUE
        )
    }

    if(gene %in% exp_gene){
        p3 <- plotEmbedding(
            ArchRProj = projHeme5, 
            colorBy = "GeneExpressionMatrix", 
            name = sort(gene), 
            embedding = "UMAP",
            quantCut = c(0.01, 0.95),
            imputeWeights = getImputeWeights(projHeme5),
            pal = paletteContinuous(set = "solarExtra")
        )
        p3_violin <- plotGroups(
            ArchRProj = projHeme5, 
            groupBy = "cell_type_plot",
            colorBy = "GeneExpressionMatrix", 
            name = sort(gene), 
            imputeWeights = getImputeWeights(projHeme5),
            plotAs="violin", size = pointSize ,
            pal = use_colors ,
            useGroups = grp_order2,
            alpha = 0.4,
            log2Norm = TRUE,
            addBoxPlot = TRUE
        )
    }

    if(gene %in% rna_gene){
        ## 表达的活性
        ## rna表达
        p4 <- FeaturePlot(scrnat, features = c(gene))

        ## rna表达的小提琴图
        p4_violin <- VlnPlot(scrnat, features = c(gene) , pt.size = 0) + FontSize(x.text = 6, y.text = 6 , x.title = 7, y.title = 7) + NoLegend()
    }

    ## 输出
    pdf(paste0(out_path, "/" , gene , "_plotEmbedding.Motif_ATAC_RNA.pdf"), width=5, height=5)
    if(gene %in% motif_gene){
        print(p1)
        print(p1_violin)
    }
    if(gene %in% atac_gene){
        print(p2)
        print(p2_violin)
    }
    if(gene %in% exp_gene){
        print(p3)
        print(p3_violin)
    }
    if(gene %in% rna_gene){
        print(p4)
        print(p4_violin)
    }
    dev.off()
}

